/*
 * accumulated_image.hpp
 *
 *  Created on: Dec 13, 2009
 *      Author: anguelos
 */

#ifndef ACCUMULATED_IMAGE_HPP_
#define ACCUMULATED_IMAGE_HPP_
#include "image.hpp"
namespace quadil{

template <typename T,typename K,int HBORDER=0,int VBORDER=0> class AccumulatedImage:public Image<T>{

public:

	AccumulatedImage(const Image<K>& in):
		Image<T>(Image<T>(in.width+HBORDER*2,in.height+VBORDER*2).getSubImage(HBORDER,VBORDER,in.width+HBORDER-1,in.height+VBORDER-1)){
		in.getAccumulativeSum(*this,HBORDER,VBORDER);

	}

	T getSum(int l,int t,int r,int b){
		return this->fp[b*this->lineStride+r]-
			   this->fp[t*this->lineStride+r]-
			   this->fp[b*this->lineStride+l]+
			   this->fp[t*this->lineStride+l];
	}
	double getAvg(int l,int t,int r,int b){
		return this->getSum(l,t,r,b)/((1.0*b-t)*(1.0*r-l));
	}


};

template <typename ACC> class Spliter{
	t_label lab;
public:
	struct Rect{
			int l,t,b,r; t_label lab;
			Rect(){
				l=0;t=0;b=0;r=0;lab=0;
			}
			Rect(int l,int t,int r,int b){
				lab=0;
				this->l=l;this->t=t;this->r=r;this->b=b;
			}
			Rect(int l,int t,int r,int b,t_label label){
				lab=label;
				this->l=l;this->t=t;this->r=r;this->b=b;
			}
			Rect(const Rect& r){

				lab=r.lab;
				this->l=r.l;this->t=r.t;this->r=r.r;this->b=r.b;
				if(r.l>r.r ||r.b<r.t){
					std::cerr<<"Rect missized: left="<<l<<" top="<<t<<" right="<<this->r<<" bottom="<<b<<"  label="<<lab<<"\n";
				}
			}
			void splitHoriz(Rect& r1,Rect& r2,int x){
				if(x>=r ||x<=l){
					std::cout<<"x="<<x<<" left="<<l<<" right="<<r<<"\n";
				}
				r1.l=l;r1.t=t;r1.r=x;r1.b=b;
				r2.l=x+1;r2.t=t;r2.r=r;	r2.b=b;
			}
			void splitVert(Rect& r1,Rect& r2,int y){
				if(y>=b ||y<=t){
					std::cout<<"y="<<y<<" top="<<t<<" bottom="<<b<<"\n";
				}
				r1.l=l;r1.t=t;r1.r=r;r1.b=y;
				r2.l=l;r2.t=y+1;r2.r=r;r2.b=b;
			}
			int getWidth(){
				return r-l+1;
			}
			int getHeight(){
				return b-t+1;
			}
			int getSurface(){
				return getWidth()*getHeight();
			}
	};
	ACC & acc;
	std::vector<Rect> leafs;
	int minW;
	int minH;
	double minScore;
	Spliter(ACC &accum):acc(accum){
		lab=1;
		minW=40;
		minH=40;
		minScore=acc.getAvg(0,0,acc.width-1,acc.height-1)/5000;
	}
	void getxSplit(Rect r,int & res,double& bestVal){
		res=(r.r+r.l)/2;
		int x;
		double avg=acc.getAvg(r.l,r.t,r.r,r.b);
		bestVal=0;
		double curVal,curLeft,curRight;
		Rect r1,r2;
		for(x=r.l+1;x<r.r-2;x++){
			r.splitHoriz(r1,r2,x);
			curLeft=acc.getAvg(r1.l,r1.t,r1.r,r1.b);
			curRight=acc.getAvg(r2.l,r2.t,r2.r,r2.b);
			curVal=((curLeft-avg)*(curLeft-avg)*(x-r.l)+(curRight-avg)*(curRight-avg)*(r.r-x))/(r.r-r.l+1);
			//curVal=(sqrt((curLeft-avg)*(curLeft-avg))+sqrt((curRight-avg)*(curRight-avg)));
			//curVal=curLeft-curRight;
			curVal=(curLeft>curRight)?curLeft:curRight;
			curVal=curVal/r.getWidth();
			if(curVal>bestVal){
				bestVal=curVal;
				res=x;
			}
		}
		res;
	}
	void getySplit(Rect r,int &res,double& bestVal){
		res=(r.b+r.t)/2;
		int y;
		double avg=acc.getAvg(r.l,r.t,r.r,r.b);
		bestVal=0;
		double curVal,curTop,curBottom;
		Rect r1,r2;
		for(y=r.t+1;y<r.b-2;y++){
			r.splitVert(r1,r2,y);
			curTop=acc.getAvg(r1.l,r1.t,r1.r,r1.b);
			curBottom=acc.getAvg(r2.l,r2.t,r2.r,r2.b);
			curVal=((curTop-avg)*(curTop-avg)*(y-r.t)+(curBottom-avg)*(curBottom-avg)*(r.b-y))/(r.b-r.t+1);
			//curVal=(sqrt((curTop-avg)*(curTop-avg))+sqrt((curBottom-avg)*(curBottom-avg)));
			//curVal=curLeft-curRight;
			curVal=(curTop>curBottom)?curTop:curBottom;
			curVal=curVal/r.getHeight();
			if(curVal>bestVal){
				bestVal=curVal;
				res=y;
			}
		}

	}
	double splitInBestDirection(Rect r,int depth){
		if(depth==0 || r.getHeight()<minH || r.getWidth()<minW){
			r.lab=lab++;
			leafs.push_back(r);
			return r.size();
		}
		Rect r1,r2;
		int y,x;
		double scorex,scorey;
		getxSplit(r,x,scorex);
		getySplit(r,y,scorey);
		if(scorex <minScore && scorey<minScore){
			return splitInBestDirection(r,0);
		}
		if(scorex>scorey){
			r.splitHoriz(r1,r2,x);
		}else{
			r.splitHoriz(r1,r2,y);
		}
		return splitInBestDirection(r1,depth-1)*splitInBestDirection(r2,depth-1);

	}
	void splitVert(Rect r,int depth){
		if(depth==0 || r.getHeight()<minH || r.getWidth()<minW){
			r.lab=lab++;
			leafs.push_back(r);
			return;
		}
		Rect r1;Rect r2;
		int y;
		double score;
		getySplit(r,y,score);
		r.splitVert(r1,r2,y);
		if(score<minScore){
			splitHoriz(r,depth-1);
			return;
		}
		splitHoriz(r1,depth-1);
		splitHoriz(r2,depth-1);
		return;
	}
	void splitHoriz(Rect r,int depth){
		if(depth==0 || r.getHeight()<minH || r.getWidth()<minW){
			r.lab=lab++;
			leafs.push_back(r);
			return;
		}
		Rect r1; Rect r2;
		int x;
		double score;
		getxSplit(r,x,score);
		r.splitHoriz(r1,r2,x);
		if(score<minScore){
			splitHoriz(r,depth-1);
			return;
		}
		splitVert(r1,depth-1);
		splitVert(r2,depth-1);
		return;
	}
	void draw(Image<t_label>& out){
		for(int i=0;i<leafs.size();i++){
			std::cout<<"Subimage("<<leafs[i].l<<","<<leafs[i].t<<","<<leafs[i].r<<","<<leafs[i].b<<") lab="<<leafs[i].lab<<"\n";
			out.getSubImage(leafs[i].l,leafs[i].t,leafs[i].r,leafs[i].b).setAll(leafs[i].lab);
		}
	}
	void print(){
		std::cout<<"Root Win=(0,0,"<<acc.width-1<<","<<acc.height-1<<")\n";
		for(int i=0;i<leafs.size();i++){
			std::cout<<"Win("<<leafs[i].l<<","<<leafs[i].t<<","<<leafs[i].r<<","<<leafs[i].b<<") lab="<<leafs[i].lab<<"\n";
		}
	}

};

}
#endif /* ACCUMULATED_IMAGE_HPP_ */
